On the melting of the nanocrystalline vortex matter in high-temperature 

superconductors 
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Multilevel Monte Carlo simulations of the vortex matter in the highly-anisotropic high- 
temperature superconductor Bi2Sr2CaCu2 0s were performed. We introduced low concentration of 
columnar defects satisfying B$ < B. Both the electromagnetic and Josephson interactions among 
pancake vortices were included. The nanocrystalline, nanoliquid and homogeneous liquid phases 
were identified in agreement with experiments. We observed the two-step melting process and also 
noted an enhancement of the structure factor just prior to the melting transition. A proposed 
theoretical model is in agreement with our findings. 

PACS numbers: 74.25.Dw, 74.25.Qt, 74.25.Ha, 74.25.Bt 



Recently there were several experiments reporting on 
some interesting features of the vortex matter melting 
transition in Bi2Sr2CaCu20s (BSCCO) in the present of 
columnar defects (CDs) [l|, 12, |3( . I n particular, these ex- 
periments consider the case that the density of the CDs 
does not exceed the density of flux-lines (FLs). Denoting 
the "matching field" by = n cc i4>o-, where n cc i is the 
number of CDs per unit area and <po = hc/2e is the flux 
quantum, we are particularly interested in < B. The 
CDs are created artificially by irradiating the sample by 
highly energetic heavy ions, like 1 GeV Pb ions, which 
produce tubes of non-conducting damaged material. The 
number of CDs can be controlled by the amount of irra- 
diation. We assume that CDs and the applied magnetic 
field are aligned with the c-axis of the single crystal sam- 
ple. The limit B ^> B^ corresponds to the formation of 
the Bose glass phase at low temperature where vortices 
are situated on the randomly placed CDs 0, 03, El Q • O n 



the other hand, for B < Ba, the picture that emerged 



from experiments 0- 12|- L| and from simulations |£J| is 
that of the "crystallites in the pores" , also referred to as 
the "porous" vortex solid. At low temperatures a skele- 
ton (or matrix) of vortices localized on CDs is formed, 
while the excess (interstitial) FLs form hexagonal crys- 
tallites in the lacuna between CDs. Thus this phase has 
a short ranged translational order, which extends to a 
distance of the order of a typical pore size. 

As the temperature is increased, and if the magnetic 
field is large enough, this heterogeneous structure melts 
in two stages: first the crystallites in the pores melt into 
a nanoliquid while the skeleton remain intact, and sub- 
sequently the skeleton melts and the liquid becomes ho- 
mogeneous. When the magnetic field is lowered (but still 
greater than B,p) these two transitions coincide into one. 
This is usually associated with a kink in the first order 
melting line. The first stage melting of the crystallites is 
observed in the experiments as a step in the equilibrium 
local magnetization 0, 0, 0] and in the simulations 
as a sharp increase in the transverse fluctuations of the 
FLs, among other signatures. The second transition was 



not observed in the experiments as a similar jump in the 
magnetization or another equilibrium property. It was 
observed in transport measurement [j| involving trans- 
port currents with alternating polarity. This establishes 
the transition as dynamical in nature. It remains to be 
established that this is also an equilibrium thermody- 
namic transition. In order to support this assertion it 
was argued Q that the second "derealization" transition 
is associated with the restoration of a broken longitudi- 
nal gauge symmetry. However, at this time it cannot be 
ruled out that the second transition is only a crossover 
associated with a gradual equilibrium change over a finite 
range of temperatures. 

The aim of this work was to try to observe the sequence 
of the two transitions and in particular the two different 
kinds of liquid phases in numerical simulations. Some 
of the advantages of using Monte Carlo simulations are 
that one can control all the interactions precisely, and 
one can take microscopic pictures and measure physical 
quantities that are difficult to measure in experiments, 
like the magnitude of the transverse fluctuations of FLs or 
the amount of their entanglement. A main disadvantage 
is that the system simulated is small and hence phase 
transitions are always broadened and are not as sharp as 
those observed in real experiments. 

The method we are using is a multilevel Monte Carlo 
simulation as reported in earlier publications Q. In a 
highly anisotropic material like BSSCO the basic de- 
grees of freedom are pancake vortices rather than string- 
like FLs. A stack of pancakes make a FL and nearest 
neighbor pancakes along the z-direction interact via the 
Josephson interaction (here we used the approximation 
recently derived in Ref. Il4l) . In addition there is a long 
range electromagnetic interaction among all pairs of pan- 
cakes that is repulsive when the two pancakes re- 
side in the same plane and attractive otherwise. Periodic 
boundary conditions are implemented in all directions. 
Similar to simulations of systems of bosons , we im- 
plement permutations when connecting FLs between the 
top and bottom planes of the simulation slab. Most sim- 
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FIG. 1: (color online) Snapshots of the FLs and CDs (green 
open circles) configurations for B = 100 G. The pancakes in 
36 planes are projected onto a single plane. Pancakes be- 
longing to the same FL are shown to be connected, (a) The 
nanocrystal phase characterized by hexagonal order in the 
pores, (b) Just before the melting. A slight enhancement of 
the structure factor (see text), (c) Nanoliquid phase. The 
pores have melted but not the skeleton, (d) Homogeneous 
liquid, full entanglement. 



ulations involved 6 x 6 x 36 pancakes making up 36 FLs in 
36 planes (some simulations involved 64 FLs). Extensive 
simulations were carried out for the case of = 50 G. 
Since we keep the number of FLs fixed at 36 this means 
that the number of CDs vary as a function of the mag- 
netic field such that N cdi = 36 B^/B. The simulations 
were carried out for fields B = 50 G to 200 G and thus 
N cd /N FL vary between 100% to 25%. We simulated be- 
tween T = 60 K to 79 K with one degree increments. 
The critical temperature for BSCCO was assumed to be 
T c = 90 K. Other parameters used were A(0) = 1700 
A (penetration depth), £(0) = 30 A (correlation length), 
d = 15 A (Cu02 plane separation), 7 = 375 (anisotropy). 
We assumed that A and £ have temperature dependence 
that is proportional to (1 - T/T^ 1 / 2 . 

An important ingredient of the simulation is determin- 
ing how to model the interaction between pancakes and 
CDs in a realistic way. There are two major sources of 
pinning: core and electromagnetic pinning. Core pinning 
|4j arises when the vortex core overlaps with a normal 
state inclusion similar to the one inside a CD. Since con- 
densation energy is lost in the vortex core, part or all 
of this energy is restored when a vortex resides inside a 
CD. Electromagnetic pinning arises ^?lll^ | when the su- 



percurrent pattern around the vortex is disturbed by the 
non-conducting defect. These two mechanisms combine 
together to yield the expression for the potential energy 
at a distance R away from the CD as felt by an individual 
pancake [lll7|: 
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where eo = <^o/(47rA) 2 is the energy scale, and r r is 
the radius of the CD. The long range tail contributing 
w —eodr 2 ,/R 2 is due mainly to the electromagnetic pin- 
ning and the short range flat region of depth w eod is 
due to the combined effects of electromagnetic and core 
pinning. If ry < V2£ a slightly different expression, 
V(R) w -0.5e dr 2 /(i? 2 + 2£ 2 ), with a similar long tail 
should be used instead 0. The importance of the long 
tail was realized recently [l3| as leading to a dependence 
B dl ~ exp(-T/T ) instead of B dl ~ exp(-T 2 /T 2 ), lead- 
ing to the decoupling of FLs from CDs at a higher tem- 
perature than previously thought . 

After experimenting with CDs of different radii 10 nm 
< r r < 30 nm we concluded that in order to obtain a good 
agreement with experimental results we needed to choose 
r r sa 30 nm in which case the formula given by Eq. Q is 
the correct one. Different heavy ions with different ener- 
gies give rise to tracks of various sizes, normally reported 
to be in the range of 4 — 20 nm in diameter. It may 
be that the actual damage caused by a track exceed its 
apparent size, thus leading to a higher effective radius. 

We now turn to the results of the actual simulations. 
Most simulations were done on the Pittsburgh super- 
computing cluster and involved averaging over 10 re- 
alizations of the disorder at each temperature. Typi- 
cally we ran 4,000,000 Monte Carlo steps for equilibra- 
tion and 2,000,000 for measurements. Snapshots are de- 
picted in Fig^at four representative temperatures. The 
nanocrystal, (a), nanoliquid (c) and homogeneous liq- 
uid (d) phases are clearly distinguishable. In addition 
we display an intermediate picture (b) taken during the 
broadened melting transition where the structure factor 
slightly increases (see later discussion of Fig. |3J| due to 
a better rearrangement of the FLs than just before the 
conclusion of the melting. In the simulations, when CDs 
are present, the transition is broadened probably by the 
fact that pores of different sizes don't melt simultane- 
ously especially in a finite size system. 

In Fig. [21 we display different measurable quantities 
that are indicators of the melting transitions. The red 
solid circles represent the pristine system with no CDs 
and the melting appears sharp both for the magnitude 
of the transverse fluctuations (a) and for the measure of 
entanglement (b), given as the fraction of FLs not end- 
ing on themselves by the periodic boundary conditions 
in the z-direction but rather, due to the proliferation of 



3 



0.4 




200 300 400 500 600 
17(1-777) (K) 




200 300 400 500 600 700 
7/(1-7/7) (K) 










-0.1 


in cake 


-0.2 


Q- 

z 


-0.3 




-0.4 




-0.5 




-0.6 




200 300 400 500 600 
7/(1-7/7) (K) 



700 



0.75 
0.5 
0.25 







0.4 










• 


• 


no CDs 






0.3 


• 




B= 150 G 






g 0.2 








♦ i \ 


[fl 

0.1 





















200 


400 


600 








T/(1-T/T 


) (K) 






















▼ B= 50 G 










♦ B = 75 G 






■ 


■ — ____ 




B = 100 G 










a B = 150G 











200 300 400 500 600 700 
7/(1-7/7) (K) 



T/(1-T/T c ) (K) 

FIG. 3: (color online) Normalized structure factor at the 
first Bragg peak for different values of the magnetic field. 
The enhancement of the structure factor just below the field- 
dependent melting transition is apparent. The inset shows 
the pristine vs. irradiated sample for B = 150 G. 



FIG. 2: (color online) Measured quantities for B — 100 G. 
The pristine case and the case of CDs with radii of 20 and 30 
nm are shown, (a) Mean square transverse fluctuations of the 
FLs. For CDs of radius 30 nm the melting and derealization 
transitions are indicated with arrows. The first stage melting 
occurs for C\ w 0.1 where Cl is the Lindemann parameter, 
(b) Number of loops of size greater than unity. A measure of 
FL entanglement, (c) Binding energy per pancake between 
FLs and CDs. (d) Number of capture pancakes per CD per 
plane. 

permutations resulting in composite loops made up of 
several individual segments. For the 30 nm CDs the two 
stage melting is quite evident. For the transverse fluc- 
tuations we first see a broad flat part ending in the first 
melting and then a second kink signaling the approach 
of the fluctuations toward the pristine value that indi- 
cates the melting of the skeleton. These transitions are 
marked by the arrows. The flat part and the final kink 
are also evident in Fig. (b) . For 20 nm the flat part is 
shorter and both transitions occur closer to the pristine 
melting. In subplots (c) and (d) we display the binding 
energy per pancake and the fraction of trapped pancakes 
per CD per plane. Both these quantities become very 
small at the derealization transition. From here onward 
we limit the discussions to CDs of radius of 30 nm. 

In Fig. |21 we display the structure factor as a func- 
tion of the reduced temperature for different fields. The 
structure factors is not a good indicator of the melting 
transition when CDs are present since the random skele- 
ton make it very low to begin with and the transition 
is broadened and becomes second order or weakly first 
order. We can see clearly the rise of the structure fac- 
tor during the onset of the melting of the nanocrystallinc 
solid since it becomes easier for the interstitial vortices 
to find better positions that are closer to hexagonal order 
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FIG. 4: (color online) Phase diagram in the T — B plane as 
suggested from the results of the simulations. 

before the full melting occurs. We found evidence for this 
by looking at snapshots at this temperature range. See 
e.g. Fig. n(b), although there are more striking exam- 
ples. A similar effect was discussed in Ref. |8| for weaker 
columnar pins of smaller radius. The inset shows the 
sharp pristine melting transition for B = 150 G as indi- 
cated by the step in the structure factor in contrast to 
the much weaker decline in the case of columnar pins. 

In Fig. 0] we display the phase diagram as it emerges 
from the results of the simulations. The green line con- 
necting the empty circles is the location of the melting 
transition for the pristine system. The blue line connect- 
ing the filled squares denote the melting of the "irradi- 
ated" systems when CDs of radius 30 nm are present. 
This line denotes the conclusion of the melting of the 
nanocrystallinc solid which now becomes a nanoliquid. 
The horizontal blue bars denote the broadening of the 
first stage melting characterized by the flat sections in 
Fig. |2 (a) and (b) . The two red lines (solid and dashed) 
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connected by horizontal error bars denote the second 
stage melting of the skeleton which can be determined 
to a precision of about 2 K. Note that for = 50 G the 
two transitions merge together. The kink in the first- 
stage melting is defined by the point of largest shift of 
the melting temperature from its pristine value and this 
occurs at B ss 100 G. A change in slope of the melting 
line is observed there. In our simulation at this point 
the two melting transitions do not yet coincide. This is 
similar to the experimental results reported in Ref. y for 
B$ — 60 G. We also show in Fig. 4 the prediction of a 
simple theoretical model described as follows: 

The cage model is often used to get a rough estimate 
of the melting transition^]. For a FL trapped by a CD 
we add the binding energy between the FL and the CD 
to the harmonic cage potential. Thus the free-energy 
functional for a single FL is given by: 



F = 



dz 



dr(z) 
dz 



U(r) 



1 



-kr 1 



(2) 



where U(r) is the binding energy 
V(r)/d as given in Eq.(l) above), e 
and k w to/a 2 , where a = y/(f> /B 
tion between FLs. We now utilize 
model to a quantum particle in a 
correspondence is h — ► T , m — > t\ 
imaginary time. In the limit L — > 
zero temperature and its energy is 
state of the corresponding Schrodin 



to the CD (U(r) = 
i is the tilt modulus 
is the mean separa- 
the mapping of this 
potential where the 
and z becomes the 
oo the particle is at 
given by the ground 
ger equation 



-^-V 2 + U(r) + hr 2 



tp(r) = Eip(r), 



(3) 



in two spatial dimensions. We have chosen the following 
value for ei 



ei « e 



1.55 + ln(A/d) 



A 2 



(4) 



described above. The derealization field obtained from 
the simulation and from the solution of the Schrodinger 
equation satisfies approximately B^i ~ exp(— T/T ) with 
To rj 3 K. This line is somewhat steeper than the experi- 
mental result , but this is also true for the melting line 
that is flatter in the experiment. This is most likely due 
to the fact that experimentally pristine samples contain 
a certain amount of point defects that lower the melting 
temperature and which are even more effective in lower 
temperatures and higher fields. 

Conclusions - In this paper we used the multilevel 
Monte Carlo method to simulate the pancake vortices 
in the highly anisotropic superconductor BSCCO in the 
presence of columnar defects of low concentration. When 
the applied field is larger than the matching field we have 
found evidence for a sequence of two melting transitions, 
one associated with the melting of crystallites and the 
second with the demise of the skeleton - the FLs attached 
to the CDs. Of course this is not a conclusive evidence 
that the latter transition is a true thermodynamic transi- 
tion. The agreement with experiment is amazingly good 
for a system of 1296 pancakes when averaging over 10 re- 
alizations of the disorder and display most of the salient 
features found in the experiments. In addition we found 
a new feature - a slight but noticeable enhancement of the 
structure factor just prior to the melting transition when 
CDs are present. This enhancement might be observed 
in experiments using small-angle neutron scattering tech- 
nique (SANS) which can effectively be used to measure 
the structure factor. 
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The first term approximatel y rep resents the contribution 
of the Josephson interaction [lj, |2(J and the second term 
results from the electromagnetic interaction ^5|. We 
solved the Schrodinger equation numerically using either 
the shooting method or alternatively the matrix method. 
The equation reduces to a one-dimensional radial equa- 
tion. The melting transition is obtained by calculating 
the expectation value of r 2 and using the Lindemann cri- 
terion with C\ rs 0.1 to identify the melting temperature. 
For U = we obtained the approximation to the pristine 
melting curve displayed in Fig. 21 Including U(r) we 
obtain the result for Bdi, the field corresponding to the 
skeleton melting. We have also solved numerically the 
Schrodinger equation without the cage and used the re- 
sulting energy and localization length to reconstruct Bdi 
following Lopatin and Vinokur 0] . The result is reason- 
able but the agreement is not as good than the method 
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